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1. Introduction 

Goodness-of-Fit (GoF) tests are designed to assess quantitatively whether a sample of 
N observations can statistically be seen as a collection of iV independent realizations 
of a given probability law, or whether two such samples are drawn from the same 
hypothetical distribution. Two well-known and broadly used tests in this class are the 
Kolmogorov-Smirnov (KS) and Cramer-von Mises (CM) tests, that both quantify how 
close an empirical cumulative distribution function (cdf) Fn is from a target theoretical 
cdf F (or from another empirical cdf) - - see Ref. [1] for a nice review and many 
references. The major strength of these tests lies in the fact that the asymptotic 
distributions of their test statistics is completely independent of the null-hypothesis 
cdf. 

It however so happens that in certain fields (physics, finance, geology, etc.) the 
random variable under scrutiny has some memory. Whereas the unconditional law of 
the variable may well be unique and independent of time, the conditional probability 
distribution of an observation following a previous observation exhibits specific patterns, 
and in particular long-memory, even when the linear correlation is short-ranged or 
trivial. Examples of such phenomena can be encountered in fluid mechanics (the velocity 
of a turbulent fluid) and finance (stock returns have small auto-correlations but exhibit 
strong volatility clustering, a form of heteroscedasticity) . The long-memory nature of 
the underlying processes makes it inappropriate to use standard GoF tests in these cases. 
Still, the determination of the unconditional distribution of returns is a classic problem 
in quantitative finance, with obvious applications to risk control, portfolio optimization 
or derivative pricing. Correspondingly, the distribution of stock returns (in particular 
the behaviour of its tails) has been the subject of numerous empirical and theoretical 
papers (see e.g. [2j [3] and for reviews [U |5] and references therein). Clearly, precise 
statements are only possible if meaningful GoF tests are available. 

As a tool to study the — possibly highly non-linear — correlations between returns, 
"copulas" have long been used in actuarial sciences and finance to describe and model 
cross-dependences of assets, often in a risk management perspective [6j [7J |5]. Although 
the widespread use of simple analytical copulas to model multivariate dependences is 
more and more criticized [HI Ej, copulas remain useful as a tool to investigate empirical 
properties of multivariate data [9]. 

More recently, copulas have also been studied in the context of auto-dependent 
univariate time series, where they find yet another application range: just as Pearson's 
p coefficient is commonly used to measure both linear cross-dependences and temporal 
correlations, copulas are well-designed to assess non-linear dependences both across 
assets or in time [TUl [HI [12] - we will speak of "self-copulas" in the latter case. 
Interestingly, when trying to extend GoF tests to dependent variables, self-copulas 
appear naturally. In our empirical study of financial self-copulas, we rely on a non- 
parametric estimation rather than imposing, for example, a Markovian structure of the 
underlying process, as in e.g. [T3| [TT]. 
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The organisation of the paper is in three parts. In Section [2] we study theoretically 
how to account for general dependence in GoF tests: we first describe the statistical 
properties of the empirical cdf of a non-iid vector of observations of finite size, as well as 
measures of its difference with an hypothesized cdf. We then study the limit properties 
of this difference and the asymptotic distributions of two norms. In Section [3] we 
go through a detailed example when the dependences are weak and described by a 
pseudo-elliptical copula. Section 0] is dedicated to an application of the theory to the 
case of financial data: after defining our data set, we perform an empirical study of 
dependences in series of stock returns, and interpret the results in terms of the "self- 
copula" ; implication of the dependences on GoF tests are illustrated for this special case 
using Monte-Carlo simulations. The concluding section summarizes the main ideas of 
the paper, and technical calculations of sections [2] and |3] are collected in the appendix. 



2. Goodness-of-fit tests for a sample of dependent draws 

2.1. Empirical cumulative distribution and its fluctuations 

Let X be a latent random vector with N identically distributed but dependent variables, 
with marginal cdf F. One realization of X consists of a time series {x\, . . . , x n , . . . , xn} 
that exhibits some sort of persistence. For a given number x in the support of F, 
let Y(x) be the random vector the components of which are the Bernoulli variables 
Y n (x) = t{ Xn <x}- The expectation value and the covariance of Y n (x) are given by: 

E\V n (x)] = F(x), (1) 
Cov(Y n {x),Y m (x')) = F nm (x,x') = C nm (F(x),F(x')), (2) 

where by definition C nm is the "copula" of the random pair (X n ,X m ). The centered 
mean of Y(x) is: 

- 1 N 

iv n=l 

which measures the difference between the empirically determined cumulative 
distribution function at point x and its true value. It is therefore the quantity on 
which any statistics for Goodness-of-Fit testing is built. Denoting u — F(x),v — F(x'), 
the covariance function of Y is easily shown to be: 

Cov (Y(u),Y(v)) = -j-(m.in(u,v) - uv) [l + V N (u,v)] (4) 



where 



q? N (u,v) = — ^ — — (5) 



LI Li 1 ! \ II. I I UV 



measures the departure from the independent case, corresponding to C nm (u, v) = uv (in 
which case ^n(u,v) = 0). Note that decorrelated but dependent variables may lead to 
a non zero value of \I/jv> since the whole pairwise copula enters the formula and not only 
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the linear correlation coefficients. When the denominator is zero, the fraction should be 
understood in a limit sense; we recall in particular that [9] 

nm{U, U) = u /i_ u } = T nm( U ) + T n. m ( l ~ U) - 1 (6) 

tends to the upper/lower tail dependence coefficients r^(l) and 7^(1) when u tends 
to 1 resp. 0. Intuitively, the presence of ty N (u,v) in the covariance of Y above leads 
to a reduction of the number of effectively independent variables, but a more precise 
statement requires some further assumptions that we detail below. 

In the following, we will restrict to the case of strong-stationary random vectors, 
for which the copula C nm only depends on the lag t = m — n, i.e. C t = C n>n+t . The 
average of A nm over n, m can be turned into an average over t: 

N-l f 

* N (u,v)=J2(l--)(A t (u,v)+A„ t (u,v)) (7) 

with A t (u,v) = A n>n+t (u, v). Note that in general A t (u,v) ^ A_ t (u,v), but clearly 
A t (u,v) = A_ t (v,u), which implies that $?n(u, v) is symmetric inu-O-u. 

We will assume in the following that the dependence encoded by A t (u,v) has a 
limited range in time, or at least that it decreases sufficiently fast for the above sum to 
converge when iV — > oo. If the characteristic time scale for this dependence is T, we 
assume in effect that T N. In the example worked out in Section [3] below, one finds: 

A t {u, v) = f (^) A< f' V \ , I(u, v) = min(w, v) - uv 

where /(•) is a certain function. If f(r) decays faster than r _1 , one finds (in the limit 
T > 1): 

. , , . , m A(u, v) + A(v, u) f°° , „, . 
(u,v) = Km * N (u,v)=T-^-f- AJ-^ / drf r , 

with corrections at least of the order of T/N when N ^> T. 



2.2. Limit properties 

We now define the process y(u) as the limit of y/~NY(u) when iV — > oo. For a given u, 
it represents the asymptotics of the difference between the empirically determined cdf 
of the underlying X's and the theoretical one, at the tt-th quantile. According to the 
Central Limit Theorem under weak dependences, it is Gaussian as long as the strong 
mixing coefficients, 

a SM (t) = supsup{|P(An B) -P(A)P(B)| : A e a({Z n (u)} n < T ), B E a({Z n (v)} n > T+t )} 

t A,B 
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associated to the sequence {Z n {u)} = {Y n {u) — u}, vanish at least as fast as 0(t~ 5 
We will assume this condition to hold in the following. For example, this condition is 
met if the function /(r) defined above decays exponentially, or if f(r > 1) = 0. 
The covariance of the process y(u) is given by: 

H(u,v)= lim N Cov (Y(u),Y(v)) = (mm(u,v) - uv) [1 + * 00 (m,u)] (8) 

JV— >oo 

and characterizes a Gaussian bridge since V[y(0)] = V[y(l)] = 0, or equivalently 
P[y(0) = y] = P[y(l) = y] = 5{y). Indeed, I(u,v) = min(-u,t>) — uv is the covariance 
function of the Brownian bridge, and ^/^(ujv) is a non-constant scaling term. 

By Mercer's theorem, the covariance H (u, v) can be decomposed on its eigenvectors 
and y(u) can correspondingly be written as an infinite sum of Gaussian variables: 



yW = E^WvV^ ( g ) 

where Zj are independent centered, unit-variance Gaussian variables, and the functions 
Uj and the numbers Xj are solutions to the eigenvalue problem: 

J H(u,v)Ui(v)dv = XiUi(u) with J Ui{u)Uj{u) du = 6 tj . (10) 

In order to measure a limit distance between distributions, a norm over the space 
of continuous bridges needs to be chosen. Typical such norms are the norm-2 (sum of 
squares, as the bridge is always integrable), and the norm-sup (as the bridge always 
reaches an extremal value). Other norms, such as the Anderson-Darling test, can also 
be considered, see Ref. pQ. 

In practice, for every given problem, the covariance function in Equation (JHJ) has 
a specific shape, since ^^(UjV) is copula-dependent. Therefore, contrarily to the case 
of independent random variables, the GoF tests will not be characterized by universal 
(problem independent) distributions. 

2.3. Law of the norm-2 (Cramer-von-Mises) 

The norm-2 of the limit process is the integral of y 2 over the whole domain: 

CM= [ yiufdu. (11a) 
Jo 

In the representation ([9]), it has a simple expression: 



CM = $>4 (116) 

| This condition means that the occurence of any two realizations of the underlying variable can 
be seen as independent for sufficiently long time between the realizations. Since the copula induces a 
measure of probability on the Borel sets, it amounts in essence to checking that \Ct (u, v) — uv\ converges 
quickly towards 0. See Refs. [HI 021 [TO] f° r definitions of a—, j3—, p~ mixing coefficients and sufficient 
conditions on copulas for geometric mixing (fast exponential decay) in the context of copula-based 
stationary Markov chains. 
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and its law is thus the law of an infinite sum of squared independent gaussian variables 
weighted by the eigenvalues of the covariance function. Diagonalizing H is thus sufficient 
to find the distribution of CM, in the form of the Fourier transform of the characteristic 
function 

<t>(t) = E [e itCM ] = J] (1 - 2it\ 3 Y^ . (12) 

3 

The hard task consists in finding the infinite spectrum of H (or some approximations, 
if necessary). 

Ordering the eigenvalues by decreasing amplitude, Equation makes explicit 



the decomposition of CM over contributions of decreasing importance so that, at a 
wanted level of precision, only the most relevant terms can be retained. In particular, 
if the top eigenvalue dominates all the others, we get the chi-square law with a single 
degree of freedom: 




P [CM < k) = erf J—. (13) 



Even if the spectrum cannot easily be determined but H(u,v) is known, all the 
moments of the distribution can be computed exactly. For example: 

E[CM] = TrH = C H(u, u) du, (14a) 
Jo 

Y[CM) =2TtH 2 = 2 f f H(u,v) 2 dudv. (Ub) 



2.4- Law of the supremum (Kolmogorov-Smirnov) 

The supremum of the difference between the empirical cdf of the sample and the target 
cdf under the null-hypothesis has been used originally by Kolmogorov and Smirnov as 
the measure of distance. The variable 

KS = sup \y(u)\ (15) 
«e[o,i] 

describes the limit behaviour of the GoF statistics. In the case where 1 + \l/ 00 (M,f) can 



be factorized as yip(u)yip(v), the procedure for obtaining the limiting distribution was 
worked out in [16], and leads to a problem of a diffusive particle in an expanding cage, 
for which some results are known. There is however no general method to obtain the 
distribution of KS for an arbitrary covariance function H. 

Nevertheless, if H has a dominant mode, the relation (Q becomes approximately: 
y(u) = U (u)y/\o~z Q = k (u )z , and 



KS=J\ \z \ sup \U (u)\ = Ko(uq)\z q \. (16) 
«e[o,i] 



The cumulative distribution function is then simply 

k 

.V2ko(«o 



F[KS <k]= erf ( k - ) , k > 0. (17) 
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This approximation is however not expected to work for small values of k, since in this 
case Zq must be small, and the subsequent modes are not negligible compared to the 
first one. A perturbative correction — working also for large k only — can be found 
when the second eigenvalue is small, or more precisely when y(u) = k,q(u)z q + k(u)z\ 
with e = k/kq <C 1. The first thing to do is find the new supremum 

M * = argsup(yH 2 )= MS + ^-^. (18) 

Notice that it is dependent upon z , z\ so that KS is no longer exactly the absolute value 
of a Gaussian. However it can be shown (after lenghty but straightforward calculation) 
that, to second order in e, y(u*) remains Gaussian, albeit with a new width 



K* ~ \J + K 2 > Kq, (19) 

where all the functions are evaluated at Uq. In fact, this approximation works also with 
more than two modes, provided 

k(u) 2 = \U 3 (u) 2 < k (u) 2 = \ U (u) 2 , (20) 

in which case: 



K * ~ pl x i u M)- (2i) 

3. An explicit example: The log-normal volatility model 

In order to illustrate the above general formalism, we focus on the specific example of 
the product random variable X = cr£, with iid Gaussian residuals £ and log-normal 
stochastic standard-deviations a = e w (again we denote generically by F the cdf of X). 
Such models are common in finance to describe stock returns, as will be discussed in 
the next section. For the time being, let us consider the case where the cu's are of zero 
mean, and covariance given by: 

Cov(u n cu n+t ) = S 2 / (|) , (t > 0). (22) 

The pairwise copulas in the covariance of Y can be explicitely written in the limit 
of weak correlations, S 2 — > 0. One finds: 

C t (u, v)-uv = S 2 / (1) A(u)A(v) + G(S 4 ) (23) 
with A(u) = / ip{u)ip\ (24) 

— oo 

where here and in the following ip(-) denotes the univariate Gaussian pdf, and $(•) the 
Gaussian cdf. The spectrum of A(u,v) = A(u)A(v) consists in a single non-degenerate 
eigenvalue \ A = TrA = Jq A(u) 2 du, and an infinitely degenerate null eigenspace. 
Assuming short-ranged memory, such that = Y,^Lif( r ) < +oo, the covariance 
kernel reads: 

H(u, v) = I(u, v) + 2TS 2 / 00 A(u, v). 
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Depending on the value of the parameters, the first term or the second term may be 
dominant. Note that one can be in the case of weak correlations (£ 2 — > 0) but long range 
memory T 3> 1, such that the product TS 2 can be large (this is the case of financial 
time series, see below). If TE 2 is small, one can use perturbation theory around the 
Brownian bridge limit (note that TrJ w lOTrA, see Appendix), whereas if TS 2 is large, 
it is rather the Brownian term I(u,v) that can be treated as a perturbation. Elements 
of the algebra necessary to set up these perturbation theories are given in the Appendix. 

It is interesting to generalize the above model to account for weak dependence 
between the residuals £ and between the residual and the volatility, without spoiling the 
log-normal structure of the model. We therefore write: 



X Q = £ e u ' ; X t = ^e^o+Aeo+VW-/^ w i t h E[£ &] = r t 
where all the variables are Af(0, 1), so that in particular 
Pt = Corr(X ,X t ) = r t {\ + /3 t 2 )e at_1 
Corr(X 2 ,X t 5 



, a v2 , (1 + 2r 2 (l + 10/3 2 + 8{3?) + 4/3 2 ) e 4 - - 1 



Corr(X ,X 2 )=2A 



3e 4 - 1 
(l + 2r 2 (l + 2/3 2 ))e 2Qt -5 



The univariate marginal distributions of Xo and X t are identical and their cdf is given 
by the integral 

/OO rv> 
<p{u)*(-)6u>. (25) 
-oo e 

Expanding the bivariate cdf (or the copula) in the small dependence parameters at, Pt, Pt 
around (0, 0, 0), we get 

C t (u, v) — uv ~ a t A{u, v) — f3 t B(u, v) + ptR(u, v) (26) 
w a t A(u)A(v) - f3 t R(u)A{v) + p t R(u)R{v) 

where A(u) was defined above in Equation §M§, and 

(p(u)tp( - { ) )du = R{l-u). (27) 

-oo e^ 1 

The contributions of A(u, v), B(u, v) and R{u, v) on the diagonal are illustrated in 
Figure [TJ Notice that the term B(u,v) (coming from cross-correlations between £ and 
uj t , i-e. the so-called leverage effect, see below) breaks the symmetry C t (u,v) ^ C t {y,u). 

We now turn to a numerical illustration of our theory, in the simple case where 
only volatility correlations are present (i.e. (3 t = pt = in Equation (|26|) above). We 
furthermore assume a Markovian dynamics for the log-volatilities: 

X n = £„e""- VM , with u n+1 = gu n + E Vn , (28) 

where g < 1 and rj n are iid Gaussian variables of zero mean and unit variance. In this 

case, 

S 2 

a t = Cov(u n u n+ t) = (29) 

1 - g* 




Figure 1. Copula diagonal of the log- normal volatility model: linear corrections to 
independence. Left: correction R(u,u) due to correlation of the residuals (vertical 
axis in multiples of p) Middle: correction A(u, u) due to correlation of the log-vols 
(vertical axis in multiples of a) Right: correction B(u,u) due to leverage effect 
(vertical axis in multiples of —ft) 



In the limit where E 2 1, the weak dependence expansion holds and one finds explicitly: 

H(u, v) = I(u, v) + 2 tf A(u, v). (30) 

In order to find the limit distribution of the test statistics, we procede by Monte-Carlo 
simulations. The range [0, l] 2 of the copula is discretized on a regular lattice of size 
(M x M). The limit process is described as a vector with M components and built 
from Equation (jUJ) as y = UA^z where the diagonal elements of A are the eigenvalues 
of H (in decreasing order), and the columns of U are the corresponding eigenvectors. 
Clearly, Cov(y,y) = UAW = H. 

For each Monte-Carlo trial, M independent random values are drawn from a 
standard Gaussian distribution and collected in z. Then y is computed using the above 
representation. This allows one to determine the two relevant statistics: 

KS = max \y u \ 

u=l...M 

i M i -i 

CM = —Yf u = — yty = _ z tAz. 

The empirical cumulative distribution functions of the statistics for a large number of 
trials are shown in Figure |2] together with the usual Kolmogorov-Smirnov and Cramer- 
von-Mises limit distributions corresponding to the case of independent variables. 

In order to check the accuracy of the obtained limit distribution, we generate 350 
series of = 2500 dates according to Equation (1251) . For each such series, we perform 
two GoF tests, namely KS and CM, and calculate the corresponding p-values. By 
construction, the p-values of a test should be uniformly distributed if the underlying 
distribution of the simulated data is indeed the same as the hypothesized distribution. 
In our case, when using the usual KS and CM distributions for independent data, the 
p-values are much too small and their histogram is statistically not compatible with the 
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Figure 2. Markovian model - Left: Cumulative distribution function of the 
supremum of y(u). Right: Cumulative distribution of the norm-2 of y(u). The 
cases of independent drawings (thin red) and dependent drawings (bold black) are 
compared. The dependent observations are drawn according to the weak-dependence 
kernel (l3"u| with parameters g — 0.88, S 2 — 0.05. Insets: The effective reduction 

ratio A / AT N , \ = ec df (") w here L = KS,CM. The dashed vertical line is located at 
V cdf L («) 

the 95-th centile and thus indicates the reduction ratio corresponding to the p-value 
p = 0.05 (as the test is unilateral). 



Kolmogorov-Smirnov Cramer-von Mises 




p-values p-values 



Figure 3. Histogram of the p-values in the GoF test on simulated data, according to 
Equation (|28|) . Uniform distribution of the p-values of a test indicates that the correct 
law of the statistics is used. 



uniform distribution. Instead, when using the appropriate limit distribution found by 
Monte- Carlo and corresponding to the correlation kernel ( I30p . the calculated p-values are 
uniformly distributed, as can be visually seen on Figure [31 and as revealed statistically 
by a KS test (on the KS test!), comparing the 350 p-values to H : p ~ U[0, 1]. 
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If, instead of the AR(1) (Markovian) prescription (|28p . the dynamics of the u n is given 
by a Fractional Gaussian Noise (i.e. the increments of a fractional Brownian motion) 
[PT] with Hurst index 2 ^ L > |, the log- volatility has a long ranged autocovariance 

at = Cov(u n , oo n+t ) = — ((t + lf- u - 2t 2 -» + \t- l\ 2 - u ) , t > (31) 

that decays as a power law oc (2 — 3v + i* 2 )t~ u as t — > oo, corresponding to long-memory, 
therefore violating the hypothesis under which the above theory is correct. Still, we 
want to illustrate that the above methodology leads to a meaningful improvement of 
the test, even in the case of long-ranged dependences. The corresponding covariance 
kernel of the Xs, 

N / t\ 

H(u, v) = I(u, v) + 2S 2 A( M , v) £ (l - W ) at, (32) 

is used in a Monte-Carlo simulation like in the previous case in order to find the 
appropriate distribution of the test statistics KS and CM (shown in Figure HJ see 
caption for the choice of parameters). We again apply the GoF tests to simulated 
series, and compute the p-values according to the theory above. As stated above, our 
theory is designed for short-ranged dependences whereas the FGN process is long-ranged. 
The p-values are therefore not expected to be uniformly distributed. Nevertheless, the 
distribution of the p-values is significantly corrected toward the uniform distribution, 
see Figure HI with the naive CM distribution (middle), the obtained p-values are 
localized around zero, suggesting that the test is strongly negative. If instead we use 
our prediction for short-range dependences (right), we find a clear improvement, as the 
p-values are more widely spread on [0, 1] (but still not uniformly). 




Figure 4. Fractional Brownian Motion - Left: Cumulative distribution function of 
the norm-2 of y(u), see Fig. [2] for full caption. Middle- Right: Histogram of the 
p-values in the CM test on simulated data, using the iid CM distribution (middle) 
and the prediction obtained assuming short range correlations (right). The dependent 
observations are drawn according to (pITj) with parameters v — |, S 2 = 1, N — 1500. 
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4. Application to financial time series 

4-1. Stylized facts of daily stock returns 

One of the contexts where long-ranged persistence is present is time series of financial 
asset returns. At the same time, the empirical determination of the distribution of these 
returns is of utmost importance, in particular for risk control and derivative pricing. As 
we will see, the volatility correlations are so long-ranged that the number of effectively 
independent observations is strongly reduced, in such a way that the GoF tests are not 
very tight, even with time series containing thousands of raw observations. 
It is well-known that stock returns exhibit dependences of different kinds: 

• at relatively high frequencies (up to a few days), returns show weak, but significant 
negative linear auto-correlations (see e.g. [18J); 

• the absolute returns show persistence over very long periods, an effect called 
multiscale volatility clustering and for which several interesting models have been 
proposed in the last ten years [T9| 1201 I2T1 [22| |23| 124] ; 

• past negative returns favor increased future volatility , an effect that goes under 
the name of "leverage correlations" in the literature [25j [26], ETJ EH [29} [30] . 

Our aim here is neither to investigate the origin of these effects and their possible 
explanations in terms of behavioral economics, nor to propose a new family of models 
to describe them. We rather want to propose a new way to characterize and measure 
the full structure of the temporal dependences of returns based on copulas, and extract 
from this knowledge the quantities needed to apply GoF tests to financial times series. 

Throughout this section, the empirical results are based on a data set consisting 
of the daily returns of the stock price of listed large cap US companies. More precisely 
we keep only the 376 names present in the S&P-500 index constantly over the five 
years period 2000-2004, corresponding to N = 1256 days. The individual series are 
standardized, but this does not change the determination of copulas, that are invariant 
under increasing and continuous transformations of the marginals. 

4-2. Empirical self-copulas 

For each (u,v) on a lattice, we determine the lag dependent "self-copula" C t (u,v) by 
assuming stationarity, i.e. that the pairwise copula C nm (u,v) only depends on the time 
lag t = m — n. We also assume that all stocks are characterized by the same self-copula, 
and therefore an average over all the stock names in the universe is done in order to 
remove noise. Both these assumptions are questionable and we give at the end of this 
section an insight on how non-stationarities as well as systematic effects of market cap, 
liquidity, tick size, etc, can be accounted for. 

The self-copulas are estimated non-parametrically with a bias correctional, then 
fitted to the parametric family of log-normal copulas introduced in the previous section. 

§ Details on the copula estimator and the bias issue are given in appendix. 
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Figure 5. Diagonal (top) and anti-diagonal (bottom) of the self-copula for different 
lags; the product copula has been subtracted. A fit with Equation (|26l) is shown in thin 
red. Note that the y scale is small, confirming that the weak dependence expansion is 
justified. The dependence is still significant even for t ~ 500 days! 



We assume (and check a posteriori) that the weak dependence expansion holds, leaving 
us with three functions of time, a t , (3 t and p t , to be determined. We fit for each t the 
copula diagonal C t (u,u) to Equation (1261) above, determine a t , (3 t and p t , and test for 
consistency on the anti-diagonal C t (u, 1 — u). Alternatively, we could determine these 
coefficients to best fit C t (u,v) in the whole (u,v) plane, but the final results are not 
very different. The results are shown in Figure [5] for lags t = 1,8,32,256 days. Fits of 
similar quality are obtained up to t = 512. 

Before discussing the time dependence of the fitted coefficients a t , fit and p t , let us 
describe how the different effects show up in the plots of the diagonal and anti-diagonal 
copulas. The contribution of the linear auto-correlation can be directly observed at the 
central point Ct(§, 5) of the copula. It is indeed known [9] that for any pseudo-elliptical 
model (including the present log-normal framework) one has: 

C tih\) = 7 + 7^arcsinp t . 
4 271 

Note that this relation holds beyond the weak dependence regime. If /3 t (B) = C t (|, |) — \ 
- this is in fact Blomqvist's beta coefficient (3TJ — the auto-correlation is measured by 
A = sin(2vrA (B) )- 

The volatility clustering effect can be visualized in terms of the diagonals of the self- 
copula; indeed, the excess (unconditional) probability of large events following previous 
large events of the same sign is (C t (u,u) — u 2 ) with u < | for negative returns, and 
u > I for positive ones. On the anti-diagonal, the excess (unconditional) probability 
of large positive events following large negative ones is, for small u < |, the upper-left 
volume (C t (u, 1) — u • 1) — (C t (u, 1—u) — u{l—u)) = u(l—u) — C t (u, 1—u) and similarly 
the excess probability of large negative events following large positive ones is the same 
expression for large u > | (lower-right volume). As illustrated on Figure these four 
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Figure 6. Auto-correlation of the volatilities, for lags ranging from 1 to 768 days. 
Each point represents the value of a t extracted from a fit of the empirical copula 
diagonal at a given lag to the relation ([26)) . We also show the fit to a multifractal 
model, at = -£ 2 log ^ with E 2 = 0.046 and T = 1467 days. 



quadrants exceed the independent case prediction, suggesting a genuine clustering of 
the amplitudes, measured by at. Finally, an asymmetry is clearly present: the effect 
of large negative events on future amplitudes is stronger than the effect of previous 
positive events. This is an evidence for the leverage effect: negative returns cause a 
large volatility, which in turn makes future events (positive or negative) to be likely 
larger. This effect is captured by the coefficient (3 t . 

The evolution of the coefficients at, 0t an d pt for different lags reveals the following 
properties: i) the linear auto-correlation p t is short-ranged (a few days), and negative; 
ii) the leverage parameter (3 t is short-ranged and, as is well known, negative, revealing 
the asymmetry discussed above; iii) the correlation of volatility is long-ranged and of 
relatively large positive amplitude (see Figure [6]), in line with the known long range 
volatility clustering. More quantitatively, we find that the parameter at for lags ranging 
from 1 to 768 days is consistent with an effective relation well fitted by the "multifractal" 
[32| |20| |33j El] prediction for the volatility autocorrelations: a t = — £ 2 log|;, with an 
amplitude X 2 = 0.046 and a horizon T = 1467 days consistent, in order of magnitude, 
with previous determinations. 

The remarkable point, already noticed in previous empirical works on multifractals 
[32| 134"] . is that the horizon T, beyond which the volatility correlations vanish, is found 
to be extremely long. In fact, the extrapolated horizon T is larger than the number 
of points of our sample N\ This long correlation time has two consequences: first, the 
parameter 2T£ 2 / 00 that appears in the kernel H(u, v) is large, ~ 135. This means 
that the dependence part TE 2 /^ A(u, v) is dominant over the independent Brownian 
bridge part I(u,v). This is illustrated in Figured where we show the first eigenvector 
of H(u,v), which we compare to the non-zero eigenmode of A(u,v), and to the first 
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u 



Figure 7. Bold black: The first eigenvector of the empirical kernel H(u, v) — 
I(u,v) [1 + ^n(u, v)}. Plain red: The function A(u) (normalized), corresponding to 
the pure effect of volatility clustering in a log-normal model, in the limit where the 
Brownian bridge contribution I(u, v) becomes negligible. Dashed blue: The largest 
eigenmode |1) = \/2sin(7nt) of the independent kernel I(u,v). 



eigenvector of I(u,v). Second, the hypothesis of a stationary process, which requires 
that iV <C T, is not met here, so we expect important preasymptotic corrections to the 
above theoretical results. 

4-3. Monte-Carlo estimation of the limit distributions 

Since H(u,v) is copula-dependent, and considering the poor analytical progress made 
about the limit distributions of KS and CM in cases other than independence, the 
asymptotic laws will be computed numerically by Monte-Carlo simulations (like in the 
example of Section [3]) with the empirically determined H(u,v). 

The empirical cumulative distribution functions of the statistics for a large number 
of trials are shown in Figure [8] together with the usual Kolmogorov-Smirnov and Cramer- 
von-Mises limit distributions corresponding to the case of independent variables. One 
sees that the statistics adapted to account for dependences are stretched to the right, 
meaning that they accept higher values of KS or CM (i.e. measures of the difference 
between the true and the empirical distributions). In other words, the outcome of a test 
based on the usual KS or CM distributions is much more likely to be negative, as it will 
consider "high" values (like 2-3) as extremely improbable, whereas a test that accounts 
for the strong dependence in the time series would still acccept the null-hypothesis for 
such values. 

As an illustration, we apply the test of Cramer- von Mises to our dataset, comparing 
the empirical univariate distributions of stock returns to a simple model of log-normal 
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Kolmogorov-Smirnov Cramer-von Mises 




Figure 8. Left: Cumulative distribution function of the supremum of y(u). Right: 
Cumulative distribution of the norm-2 of y(u). The cases of independent drawings (thin 
red) and dependent drawings (bold black) are compared. The dependent observations 
are drawn according to the empirical average self-copula of US stock returns in 2000- 

2004. Insets: The effective reduction ratio ,/ .. N , , = ec . d L (u) where L = KS,CM. 

V N eSW cdfy («) 



stochastic volatility 

X = e s "- s2 £ where ^,w~JV(0,l). (33) 
The volatility of volatility parameter s can be calibrated from the time series {x t }t as 

We want to test the hypothesis that the log-normal model with a unique value of s for 
all stocks is compatible with the data. In practice, for each stock i, is chosen as the 
average of (I34p over all other stocks in order to avoid endogeneity issues and biases in 
the calculations of the p-values of the test. s, is found to be ~ 0.5 and indeed almost 
identical for all stocks. Then the GoF statistic CM is computed for each stock i and 
the corresponding p-value is calculated. 

Figure M shows the distribution of the p-values, as obtained by using the usual 
asymptotic Cramer-von Mises distribution for independent samples (left) and the 
modified version allowing for dependence (right). We clearly observe that the standard 
Cramer-von Mises test strongly rejects the hypothesis of a common log-normal model, as 
the corresponding p-values are concentrated around zero, which leads to an excessively 
high rejection rate. The use of the generalized Cramer-von Mises test for dependent 
variables greatly improves the situation, with in fact now too many high values of p. 
Therefore, the hypothesis of a common log-normal model for all stocks cannot be rejected 
when the long- memory of volatility is taken into account. The overabundant large values 
of p may be due to the fact that all stocks are in fact exposed to a common volatility 
factor (the "market mode"), which makes the estimation of s somewhat endogeneous 
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Figure 9. Histogram of the p-values in a Cramer-von Mises-like test, see the text for 
the test design. Left: when using the standard Cramer-von Mises law, the obtained 
p-values are far from uniformly distributed and strongly localized under the threshold 
p = 0.05 (dashed vertical line) occasioning numerous spurious rejections. Right: 
when using the modified law taking dependences into account, the test rejects the 
hypothesis of an identical distribution for all stocks much less often. 



and generates an important bias. Another reason is that the hypothesis that the size 
of the sample N is much larger than the correlation time T does not hold for our 
sample, and corrections to our theoretical results are expected in that casejjj] It would 
be actually quite interesting to extend the above formalism to the long-memory case, 
where T > N > 1. 

4-4- Beyond stationarity and universality 

The assumption that financial time series are stationary is far from granted. In fact, 
several observations suggest that financial markets operate at the frontier between 
stationarity and non-stationarity. The multifractal model, for example, assumes 
that correlations are decaying marginally slowly (as a logarithm), which technically 
corresponds to the boundary between stationary models (when decay is faster) and 
non-stationary models. Furthermore, as stated above, the horizon T that appears in the 
model is empirically found to be very long (on the order of several years) and therefore 
very difficult to measure precisely. Other models, like the "multi-scale" GARCH |22j EH] , 
point in the same direction: when these models are calibrated on financial data, the 
parameters are always found to be very close to the limit of stability of the model (see 
|35j for a discussion of this point). 

Furthermore, what is relevant in the present context is strong stationarity, i.e. the 
stationarity of the full self-copula. Hence, testing for stationarity amounts to comparing 
copulas, which amounts to a GoF for two-dimensional variables . . . for which general 
statistical tools are missing, even in the absence of strong dependence! So in a sense 

|| Note that in practice, we have estimated ^n(u, v) by summing the empirically determined copulas 
up to £ max = 512, which clearly underestimates the contribution of large lags. 
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this problem is like a snake chasing its own tail. A simple (and weak) argument is 
to perform the above study in different periods. When fitting the multifractal model 
to the self-copulas, we find the following values for (E 2 ,lnT): (0.032,6.39) in 1995- 
1999, (0.046,7.29) in 2000-2004, (0.066,8.12) in 2000-2009 and (0.078,7.86) in 2005- 
2009. These numbers vary quite a bit, but this is exactly what is expected with the 
multifractal model itself when the size of the time series N is not much larger than the 
horizon T! (see 136}|37] for a detailed discussion of the finite N properties of the model). 

As of the universality of the self-copula across the stocks, it is indeed a reasonable 
assumption that allows one to average over the stock ensemble. We compared the 
averaged self-copula on two subsets of stocks obtained by randomly dividing the 
ensemble, and found the resulting copulas to be hardly distinguishable. This can be 
understood if we consider that all stocks are, to a first approximation, driven by a 
common force — the "market" — and subject to idiosyncratic fluctuations. We know 
that this picture is oversimplified, and that systematic effects of sector, market cap, 
etc. are expected and could in fact be treated separately inside the framework that we 
propose. Such issues are related to the cross-sectional non-linear dependence structure 
of stocks, and are far beyond the objectives of this article. 

5. Conclusion 

The objectives of this paper were twofold: on the theoretical side, we introduced a 
framework for the study of statistical tests of Goodness-of-Fit with dependent samples; 
on the empirical side, we presented new measurements as well as phenomenological 
models for non-linear dependences in time series of daily stock returns. Both parts 
heavily rely on the notion of bivariate self-copulas. 

In summary, GoF testing on persistent series cannot be universal as is the case for 
iid variables, but requires a careful estimation of the self-copula at all lags. Correct 
asymptotic laws for the Kolmogorov-Smirnov and Cramer-von Mises statistics can be 
found as long as dependences are short ranged, i.e. T <C N. From the empirical 
estimation of the self-copula of US stock returns, long-ranged volatility clustering with 
multifractal properties is observed as the dominant contribution to self-dependence, in 
line with previous studies. However, subdominant modes are present as well and a 
precise understanding of those involves an in-depth study of the spectral properties of 
the correlation kernel H. 

One of the remarkable consequences of the long-memory nature of the volatility 
is that the number of effectively independent observations is significantly reduced, as 
both the Kolmogorov-Smirnov and Cramer-von Mises tests accept much larger values of 
the deviations (see Figure EJ). As a consequence, it is much more difficult to reject the 
adequation between any reasonable statistical model and empirical data. We suspect 
that many GoF tests used in the literature to test models of financial returns are 
fundamentally flawed because of the long-ranged volatility correlations. In intuitive 
terms, the long-memory nature of the volatility can be thought of as a sequence of 
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volatility regime shifts, each with a different lifetime, and with a broad distribution of 
these lifetimes. It is clear that in order to fully sample the unconditional distribution of 
returns, all the regimes must be encountered several times. In the presence of volatility 
persistence, therefore, the GoF tests are much less stringent, because there is always a 
possibility that one of these regimes was not, or only partially, sampled. The multifractal 
model is an explicit case where this scenario occurs. 

We conclude with two remarks of methodological interest. 

1) The method presented for dealing with self-dependences while using statistical 
tests of Goodness-of-Fit is computationally intensive in the sense that it requires to 
estimate empirically the self-copula for all lags over the entire unit square. In the non- 
parametric setup, discretization of the space must be chosen so as to provide a good 
approximation of the continuous distance measures while at the same time not cause 
too heavy computations. Considering that fact, it is often more appropriate to use the 
Cramer- von Mises-like test rather than the Kolmogorov-Smirnov-like, as numerical error 
on the evaluation of the integral will typically be much smaller than on the evaluation 
of the supremum on a grid, more so when the grid size is only about -g « 

2) The case with long-ranged dependence T ^> N ^> 1 cannot be treated in the 
framework presented here. First because the Central Limit Theorem does not hold 
in that case, and finding the limit law of the statistics may require more advanced 
mathematics. But even pre-asymptotically, summing the lags over the available data 
up to t w iV means that a lot of noise is included in the determination of *&n(u,v) 
(see Eqation E]). This, in turn, is likely to cause the empirically determined kernel 
H(u,v) not to be positive definite. One way of addressing this issue is to follow a 
semi-parametric procedure: the copula Ct is still estimated non-parametrically, but the 
kernel H sums the lagged copulas C t only up to a scale where the linear correlations 
and leverage correlations vanish, and only one long-ranged dependence mode remains. 
This last contribution can be fitted by an analytical form, that can then be summed up 
to its own scale, or even to infinity. 

In terms of financial developments, we believe that an empirical exploration of the 
self-copulas for series of diverse asset returns and at diverse frequencies is of primordial 
importance in order to grasp the complexity of the non-linear time-dependences. In 
particular, expanding the concept of the self-copula to pairs of assets is likely to 
reveal subtle dependence patterns. From a practitioner's point of view, a multivariate 
generalization of the self-copula could lead to important progresses on such issues as 
causality, lead-lag effects and the accuracy of multivariate prediction. 
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Appendix A. Pseudo-elliptical copula: expansion around independence 



We compute here the spectrum and eigenvectors of the kernel H(u, v) in the case of 
pseudo-elliptical copula with weak dependences, starting from the expansion (126]) . 

The situation is better understood in terms of operators acting in the Hilbert space 
of continuous functions on [0, 1] vanishing in the border. Using Dirac's braket notations, 
A = \A){A\, B = \R)(A\, R = \R)(R\. The sine functions \j) = y/2 sin(jvra) build a 
basis of this Hilbert space, and interestingly they are the eigenvectors of the independent 
kernel I(u, v) (I stands for 'Independence' and is the covariance matrix of the Brownian 
motion: I = M — P where M denotes the bivariate upper Frechet-Hoeffding copula and 
P the bivariate product copula). 

It is then easy to find the spectra: rank-one operators have at most one non- 
null eigenvalue. Using the parities of A(u) and R(u) with respect to ~ and imposing 
orthonormality of the eigenvectors, we can sketch the following table of the non zero 
eigenvalues and eigenvectors of the different operators: 

Aj = (jvr)- 2 Uj(u) = \j) 

\ R =(R\R) = TrR \U*) = \R)/y/TrR 
X A =(A\A) = TrA \U£) = \A)/VT^A 

For the pseudo-elliptical copula with weak dependence, H has the following general 
form: 

H = I + pR + aA-^(B + B^). (A.l) 

The operator B + has two non zero eigenvalues ±a/ \ R \ A , with eigenvectors 
(\Uq) ±\Uq))/ a/2. In order to approximately diagonalize H, it is useful to notice that in 
the present context A and R are close to commuting with /. More precisely, it turns out 
that \Uq) is very close to |2), and \Uq) even closer to Indeed, a 2 = (U^\2) w 0.9934 
and T\ = (Uq"\1) ~ 0.9998. Using the symmetry of A and R, we can therefore write: 

|[/ A ) = a 2 |2)+e a |2 ± ) with (2|2 ± > = (2j-l|2 ± ) = 0, Vj > 1 
\U*) =n|l) + e r |l ± ) with (1|1 ± ) = (2j\l ± ) = 0,Vj > 1 

where e a = Jl — a| <^ 1 and e r = \Jl — r\ <^ 1. The components of \2±) on the even 
eigenvectors of I are determined as: 

(W)^ i>2, 

Table Al. Traces of the operators appearing in the covariance functions (multiples 
of 10~ 2 ). Traces of the powers of the rank-one A, R equal powers of their traces. The 
trace of B + B< is zero. 
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and similarly: 

(1±|2J _ 1} = mizA ,> 2 . 

Using the definition of the coefficients a t , (3 t and p t given in section 3, we introduce 
the following notations: 

N-l 

at 



/ t \ 
a = 2 TrA lim V 1 

JV^oo N J 



JV_1 / t\ 

p = 2 TiR lim V 1 - — p t 

/ N-l , , s 

/3 = 2 VTrATri? lim V 1 & 

so that H writes: 

H = I +a\U^)(U^\j-p\U^)(U^\-Mui){U^\ 
= H + e a (aa.2 |2H2 ± | -^ lfli |1)(2 ± |) 

+ e r (pn |lIT(l| -^a 2 |lIT(2|) 



+ (ae 2 a \2 ± )(2 ± \+ pe 2 r \l ± )(l ± \ - pe a e r |2j_7(i±|) 



where IV'iKV^H ^OV'iKV^I + iV^KV'il) an d -^o is the unperturbed operator (0-th order 
in both es) 

#o = E + + «^)|2)(2| + (A{ + P >?)|1)(1| - Pna 2 |2K1| 

i>3 

the spectrum of which is easy to determine as: 

A fo= A _ A / |tffo )= __H = ^4° |l) 

v H-) 



(+1+) 

Af=Aj |C/f }=|j) (j > 3) 

where 

A{ + ~pr\ + A 2 + aa\ ± y^Af + ~pr\ -\{- aa 2 2 f + 4(/3na 2 ) 2 

A± = 2 

and |±) the corresponding eigenvectors, which are linear combination of |1) and |2) 
only. Therefore, (lj_|±) = (2j_|±) = 0. This implies that there is no corrections to the 
eigenvalues of H to first order in the es. 

At the next order, instead, some corrections appear. We call: 

Vij = (priW? ) - ^(2\U^))(j\l ± )e r 
+ {aa 2 (2\U?°)-^(l\U?°))(j\2 A _)e a 
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the matrix elements of the first order perturbation of H, whence 
Af = Af + , Ho 1j x h 

j>3 A l _ A j 

V 2 

A 2 - A 2 "I" X H _ y H 

j>3 A 2 A j 

Af = Af + E . g0 ^. g „ + (&£(j|2±) a + P^(jIU) 2 - Pe a e r (j\W(j\W) 

As of the eigenvectors, it is enough to go to first order in es to get a non-trivial 
perturbative correction: 

j>3 A l - A j 
j>3 A 2 - A j 

wf) = w + £ tAjIO 

i=l,2 A j A i 

The special case treated numerically in section 3 corresponds to p = {3 = 0, such that 
the above expressions simplify considerably, since in that case Vij = and V^j-i = 0, 
while V2,2j = & a 2(Uo'\2j) . To first order in the es, the spectrum is not perturbed and 
calls Af = Af° = Af + QOjfe, so that the characteristic function of the modified CM 
distribution is, according to Equation (112"]) . 



= n (i - 2it/(j 



7T X 



1 - 2itA^ 
^ 1 - 2it Af ' 



Its pdf is thus the convolution of the Fourier transform of <f>i(t) (characteristic function 
associated to the usual CM distribution [16]) and the Fourier transform of the correction 



4> c {t) = yl — 2iiA|/y 1 — 2i£Af°. Noting that (1— 2ia 2 t)~2 is the characteristic function 
of the chi-2 distribution, it can be shown that for k > 0, and with ji = Af° for the sake 
of readability: 

' FT (0 C ) = 5(k) - [" dA-^ ( X 2 (A:; ft) * X 2 (k; A)) 



^tt v 7 J\{ dk 

V+u 



- /"dA^-^ ((/.- A)Ji(^fc) - 0* + A)Jo(^A 
8(A/i)5 V 4 V 4A At 



__fc_ aa\ aa 2 



where x 2 (k',o~ 2 ) = (2na 2 ke k ^ a2 )^^ is the pdf of the chi-2 distribution, I n are the 
modified Bessel functions of the first kind, and * denotes the convolution operation. 
The approximation on the last line holds as long as a <C X 2 = {2-k)~ 2 and in this regime 
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we obtain finally 

V[CM = k} = v / 27FT(0)(A;) = (FT(0j) * FT(0 c ))(Jfe) 

= P 7 (A;) + 4aaV /* VAz)^^ UAaal^Uk - z))Az 
Jo 

= VAk) + Aaajn 4 [ VAk - z)e- 2 ^ z I Q (AaaWz)dz 
J o 
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Figure Al. Bias of the non-parametric naive estimator of the independence 
copula with unknown marginals, for different sample sizes N = 251 (blue), 1256 
(red), 6280 (black). Top: Relative bias nv ~~ ^ Bottom: Absolute bias 

(t^TT^T - X ) /(min(w, v) - uv) 



uv 
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Appendix B. Non-parametric copula estimator 

The copula C(u,v) of a random pair (X, Y) is 

C(u,v) = P[F X (X) < u,F Y (Y) <v} = E [t {x < F ^ {u)} t {Y < F ^ {v)} 

If the univariate marginals F X ,F Y are known, the usual empirical counterpart to the 
expectation operator can be used to define the empirical copula over a sample of N i.d. 
realisations (JQ, Yj) of the random pair: 

_ 1 N 

C(U,V) = J7T, 1 {X l <F x 1 (u)} 1 {Y l <F- 1 (v)} 
iv i=l 

which is clearly unbiased for any N. But if the marginals are unknown, they have 
to be themselves estimated, for example by their usual empirical counterpart. Since 
F(x) = ~P[X < x] = E[l{x<a;}], an unbiased ecdf is obtained as F(x) = J2iLi ^{x,<x}, 
but the expectation value of 

N 

JI {F x (x l )<«} J1 {Fy(y l )<i,} 



C 'K V ) = ^E 1 {F X (X 1 )<«} 1 {F, 
i=l 



is not C{u,v) anymore (only asymptotically is C(u,v) unbiased), but rather 

E\C(u,v)\ = J dF XY (x,y)F[B(x) < Nu-l,B(y) < Nv - 1] 

where B x (x) = J2j<N ^-{x^x} has a binomial distribution B(p, N — l) with p = F(x) and 
is not independent of B Y {y). As an example, the expected value of the independence 
(product) copula C(u,v) = uv is 

resulting in a relative bias n(u,v) — 1 = ■^^■^p — 1 vanishing only asymptotically. 
As E[C(-u,f)] may not be computable in the general case, we define as 
Nu Nv 

[Nu\ [Nv\ ° M 

our non-parametrical estimator of the copula with bias correction, even when C is not 
the independence copula. Therefore, our estimator is technically biased at finite N but 
with a good bias correction, and asymptotically unbiased. 
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